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Abstract. Random networks with specified degree distributions have been pro- 
posed as realistic models ol population structure, yet the problem ol dynamically 
modeling SIR-type epidemics in random networks remains complex. I resolve this 
dilemma by showing how the SIR dynamics can be modeled with a system of three 
nonlinear ODE's. The method makes use of the probability generating function 
(PGF) formalism for representing the degree distribution of a random network 
and makes use of network-centric quantities such as the number of edges in a 
well-defined category rather than node-centric quantities such as the number of 
infecteds or susceptibles. The PGF provides a simple means of translating be- 
tween network and node-centric variables and determining the epidemic incidence 
at any time. The theory also provides a simple means of tracking the evolution of 
the degree distribution among susceptibles or infecteds. The equations are used 
to demonstrate the dramatic effects that the degree distribution plays on the fi- 
nal size of an epidemic as well as the speed with which it spreads through the 
population. Power law degree distributions are observed to generate an almost 
immediate expansion phase yet have a smaller final size compared to homoge- 
neous degree distributions such as the Poisson. The equations are compared to 
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stochastic simulations, which show good agreement with the theory. Finally, the 
dynamic equations provide an alternative way of determining the epidemic thresh- 
old where large-scale epidemics are expected to occur, and below which epidemic 
behavior is limited to finite-sized outbreaks. 

1. Introduction 

Contact patterns constitute an important aspect of heterogeneity within a popu- 
lation of susceptible and infectious individuals, but it has been a difficult factor to 
incorporate into epidemiological models. Compartment models can capture many 
aspects of population heterogeneity, such as with respect to heterogeneous suscep- 
tibility and infectiousness 26 1,6, however such models usually assume individ- 
uals mix homogeneously within each category. In contrast, the contact patterns 
responsible for the spread of many infectious diseases tend to be characterized by 
constant relationships marked by high levels of heterogeneity in the number of 
contacts per individual. 

An alternative approach is to model a population of susceptibles and infecteds 
and the contact patterns among them as a static random network 13,25,20, 
[9]. This approach has generated a new category of epidemiological models in 
which epidemics spread from node to node by traversing network connections [231 
14, 19, 27,5,24 . Random networks with specified degree distributions have been 
proposed as a simple but realistic models of population structure. This case has 
the advantage of being well understood mathematically. The expected final size 
of epidemics in random networks with a given degree distribution has been solved 
exactly [141119] . The network approach has the advantage that the mathematics of 
stochastic branching processes 28,11,2 can be brought to bear on the problem. 
This allows for precise descriptions of the distribution of outbreak sizes early in 
the course of the epidemic as well as the final size. 14, 19 
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A shortcoming of the network model is that it is difficult to describe the explicit 
dynamical behavior of epidemics on networks. The distribution of outbreak sizes 
is easy to calculate, yet the dynamic epidemic incidence, that is the number of 
infecteds at a time t, has been difficult to derive. Simulation has been used in this 
case [8]. 

Heterogeneity in the number of contacts within networks makes it difficult 
to derive differential equations to describe the course of an epidemic. Neverthe- 
less, several researchers [3, 21 , 22,4, 7] have been successful modeling many of the 
dynamical aspects of network epidemics, particularly in the early stage where 
asymptotically correct equations for disease incidence are known. These solutions 
break down, however, when the finite size of a population becomes a significant 
factor. We improve upon these results by presenting a system of nonlinear ordi- 
nary differential equations which can be used to solve for epidemic incidence at 
any time, from an initial infected to the final size, as well as other quantities of 
interest. We treat the simplest possible case of the SIR dynamics with constant 
rate of infection and recovery. Section [2] describes the model. Several examples are 
given in section [3] and section |3.1| compares the analytical results to stochastic 
simulations. 

2. SIR in Random Networks 

The networks considered here are random networks with an arbitrary degree dis- 
tribution pk (pk being the probability of a random node having degree k) [20II18| . 
Nodes can be in any of three exclusive states: susceptible (S), infectious (I), or 
recovered (TV). The dynamics are as follows. When a node is infectious, it will 
transmit infection to each of its neighbors independently at a constant rate r. 
Infectious nodes become recovered at a constant rate fj., whereupon they will no 
longer infect any neighbors. This will be made precise in the next section. 
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It is desirable to determine the dynamics of the number of susceptibles and 
infecteds and to develop equations in terms of those quantities. This, however, 
turns out to be intractable due to heterogeneity in the number of contacts. The 
problem can be resolved by developing equations in terms of dynamic variables 
representing network-based quantities, for example, the number of connections to 
susceptible or infectious nodes at a time t. The network- and node-based quantities 
are defined in the next section. 

To bridge the divide between connection- and node-based quantities, a mathe- 
matical device known as a probability generating function (PGF) |28| is extremely 
useful. The PGF has many useful properties and is frequently used in probability 
theory and the theory of stochastic branching processes. Given a discrete proba- 
bility density pt, the PGF is defined as the series: 



The variable x in the generating function serves only as a place-holder. To il- 
lustrate the utility of this device, consider the possibility that the probability of 
a node being infected, say A, is compounded geometrically according the node's 
degree. Then, the probability of a degree k node being susceptible is (1 — A) fe , that 
is, the probability of not being infected along any of k connections. If the hazard 
is identical for all nodes, the cumulative epidemic incidence (the fraction of nodes 
infectious or recovered) will be 



g(x) = pa + p±x + p 2 x 2 + psx 3 H 



(1) 



J = 1 — [po(l — A) +pi(l - A) 1 + p 2 (l - A) 2 + ■ ■ ■ ] 



(2) 



1-3(1 -A) 



(3) 



Table [T] gives a summary of the parameters used in the model. 



2.1. Definitions 

An undirected network can be defined as a graph Q = {V, £} consisting of a set of 
vertices V corresponding to the nodes in the network, and a set of edges £ with 
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Table 1. Parameters and dynamic variables for the network SIR model. 

— r :— Force of infection. The constant rate at which infectious nodes infect a 
neighbor. 

— fi := Recovery rate. The constant rate at which infected nodes become recov- 
ered. 

— pk := The probability that a node will have degree k. 

— g(x) :— The probability generating function for the degree distribution pt- 

— S := The fraction of nodes susceptible at time t. 

— I := The fraction of nodes infectious at time t. 

— R :— The fraction of nodes recovered at time t. 

— J — I + R The cumulative epidemic incidence at time t. 

— Ax Set of arcs (ego, alter) such that node ego is in set X. 

— Mx Fraction of arcs in set Ax ■ 

— Axy Set of arcs (ego, alter) s.t. ego £ X and alter G Y. 

— Mxy Fraction of arcs in set Axy- 

elements of unordered pairs of vertices, {a, b} where a, b G V. Two vertices a, b 
are said to be neighbors or neighboring each other or simply connected if there 
exists an edge e = {a, b} G £ . For the purposes of this model, the terms "vertex" 
and "node" will often be used interchangeably. 

For the random networks considered here, the probability of being connected 
to a node is proportional to the degree of that node. Denote the degree of a 
node v G V as dy Then given an edge {ft, x\ G S, the probability that x — b 
is db/ ~}2 iGV di. This definition allows multiple edges to the same node as well as 
loops from a node to itself, however the existence of multiple edges and loops is 
exceedingly rare for large sparse random networks such that results based on this 
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case can be safely applied to networks without multiple edges. Networks of this 
type can be generated by a variatiorQ of the "configuration model" [17] : 

1. To each node v € V assign an i.i.d. degree 8 V from distribution 

2. Generate a new set X of "half-edges" with S v copies of node v for all nodes 

3. Insure X has an even number of elements, for example, by deleting a uniform 
random element if odd. 

4. While X is not empty, draw two elements Vi , i>2 uniformly at random and 
create edge {«i,«2}- 

At any point in time, a vertex can be classified as susceptible, infectious, 
or recovered. Let S,T, and 7Z denote the disjoint sets of vertices classified as 
susceptible, infectious, or recovered respectively. J — I U 7Z will denote the set 
of infectious or recovered nodes. S, I, and R will denote the fraction of nodes in 
the sets S,I, and 1Z respectively. The cumulative epidemic incidence will be the 
fraction of nodes in set J . 

As stated in the previous section, infectious vertices a £ I will infect neighbor- 
ing susceptible vertices b £ S at a constant rate r. Infectious vertices will become 
recovered (move to set 1Z) at a constant rate fi. 

Although the network is undirected in the sense that any two neighboring 
vertices can transmit infection to one another, we wish to keep track of who 
infects who. Therefore, for each edge {a, b} € £ , let there be two arcs, which will 
be defined to be the ordered pairs (o, b) and (b, a). Let A denote the set of all 
arcs in the network. The first element in the ordered pair (a, b) will frequently be 
called the ego and the second element the alter. 

Axy will denote the subset of arcs such that ego € X and alter G Y. Ax will 
denote the subset of arcs such that ego £ X. Mxy = #{Axy} / '#{.4} will denote 
the fraction of arcs in the corresponding set Axy- 

1 Note that this version of the configuration model allows loops and multiple- 
edges. 
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For example, two variables will be especially important in the derivations that 
follow. Mss is the fraction of arcs with a susceptible ego and a susceptible alter. 
Msi is the fraction of arcs with a susceptible ego and and infectious alter. Ms 
will be the fraction of arcs with a susceptible ego and an alter of any type. 

2.2. Dynamics 

Our objective is to develop a deterministic model to describe epidemic dynamics 
expressed with a low-dimensional system of differential equations. At first, this 
goal may seem incompatible with network-SIR dynamics described in the last 
section. Infection spreads along links in a random network, which implies the 
epidemic incidence at any time as well as the final size must also be random, 
depending on the particular structure of a given random network. This is true, 
however it is possible to avoid such considerations by focusing on epidemic dy- 
namics in the limit as population size goes to infinity. This strategy has been used 
in previous work to calculate the expected final size of epidemics in infinite ran- 
dom networks [19] expressed as a fraction of the total population size. A similar 
strategy is followed here by considering the fraction of nodes in sets S,J, and 1Z, 
after a small fraction e nodes are infected initially in a susceptible population. 
The conclusion is the system of equations given in table [3] in terms of the dy- 
namic variables given in table [2] The dynamics predicted by these equations are 
compared to stochastic simulations with large but finite networks in section [37T| 
Consider a susceptible node ego at time t with a degree k. Then there will be 
a set of k arcs {{ego, alter\), (ego, alter 2), ■ ■ ■ , {ego, alter t)} corresponding to ego. 
We will assume that for each arc (ego, alter i) there will be a uniform probability 
pi = Msi /Ms that alteri is infectious. Then there is an expected fraction kpi arcs 
(ego, alter) such that alter is infectious. In a time dt, an expected number rkpi dt 
of these will be such that the infectious alter transmits to ego. Consequently, the 
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Table 2. Network-based dynamic variables for the network SIR model. 

— 9 :— The fraction of degree one nodes that remain susceptible at time t. 

— pi :— Ms i /Ms- The probability that an arc with a susceptible ego has an 
infectious alter. 

— ps '■= Mss/Ms- The probability that an arc with a susceptible ego has a 
susceptible alter. 

hazard for ego becoming infected at time t is 

\ k (t)=rk Pl (t) (4) 

Now let Uk (t) represent the fraction of degree k nodes that remain susceptible 
at time t, or equivalently the probability that ego in the previous example is 
susceptible. Using equation [4j 

fifc(t) = exp{- / A fc (r)dr} = exp{- / rkpi(r)dT} 

= exp{— / rp/(r)dr} fc 



Subsequently we will use the symbol 9 to denote m = exp{— f rpi(r)dT}. 
From equation [H] it is clear that Uk = 9 k . 

Given 9, it is easy to determine the fraction of nodes which remain susceptible 
at a time t. 

S = po + pitti + P2U2 + P3U3 ■ ■ ■ 

(6) 

= Pi9 + p 2 9 2 + p a 9 a + ■ ■ ■ =g(9) 

This equation makes use of the generating function g(-) for the degree distribution 
which greatly simplifies this and subsequent equations. 
The dynamics of 9 are dependent on the hazard Ai. 

d9/dt , , , 

9 (7) 
9 = -0A(t) = -6 r PI 
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Unfortunately, this does not completely specify the dynamics of 9 and by 
extension S, which also depends on the variable pi . The derivation of the dynamics 
of pi follows. 

. _ d_ Msi_ _ Msi_ _ M s M s i , , 

PI ~ dt M s ~ Ms M 2 S 1 J 

Our goal is to put equation [8] in terms of the variables 6,ps,pi and the PGF 
g(-). Ms is easily placed in terms of these variables. 

Ms = ^^Pfe x k x Pr [degree k node susceptible] / kpk 



k 

' d 



k 

Msi follows easily. 



dJ^ 



(9) 



lg\l)=eg'(e)/g\l) 

x=l 



Msi = M s x M S i/M s = M s pi = piOg (0)/g'(l) (10) 
In time dt, —S nodes become infectious. Since S = g(9), 

s = J t s= Jt m = ° 9>{e) = -^^w ( u ) 

Calculating Msi requires careful consideration of the rearrangement of arcs 
among sets Ass and Asi as — S nodes become infected in a small time interval. 
Since the hazard of becoming infected is proportional to the number of arcs to an 
infectious alter, a newly infected node will be selected with probability propor- 
tional to the number of arcs from the node to infectious nodes. 

To clarify subsequent calculations, I will introduce the notation Sxy to repre- 
sent the average degree of nodes in set X, selected with probability proportional 
to the number of arcs to nodes in set Y, not counting one arc to nodes of type 
Y. For example, if we select an arc {ego G X, alter G Y) uniformly at random 
out of the set of arcs from nodes in set X to nodes in set Y (Axy), and follow it 
to the node in set X, (ego), then Sxy will represent the average number of arcs 
(ego, alter') not counting the arc we followed to ego. This is commonly called the 
"excess degree" of a node [15]. Furthermore, Sxy(Z) will be as Sxy but counting 
only arcs from ego to nodes in set Z, (ego, alter G Z). 
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To calculate Msi we need to first calculate 5s i, and for this it is necessary to 
derive the degree distribution among susceptible nodes. It is necessary to assume^] 
that arcs from a susceptible ego to nodes in sets S,X,1Z are distributed multino- 
mially with probabilities ps,pi, and p R — 1 — ps — pi respectively. Let d sgo (X) 
be the r.v. denoting the number of arcs from ego to nodes in set X. Letting c 
normalize the distribution, and letting the dummy variables xs,Xi, and x R cor- 
respond to the number of arcs from a susceptible ego to an alter in sets S,I,1Z 
respectively, the degree distribution for susceptible nodes will be generated by 
gs(xs,xi,x R ) = ^2p k u k xsx^x^'^x^S) = i,d(I) = j\ps,Pi]/c 

(12) 

Using the multinomial theorem this becomes 

gs(xs,xi,x R ) = ^2pk0 k (x s ps +xipi +x R (l -p s -pi)) k /c 

k (13) 

= g(0(x s ps + xipi + xr(1 - ps - pi)))/g{9), 
where c = ~^Z k Pk8 k (ps + Pi + (1 — PS — Pi)) k = <?(#) normalizes the distribution. 

The degree distribution for susceptible nodes selected with probability propor- 
tional to the number of arcs to infectious nodes will be generated by the following 
equation. Note that this equation does not count one arc to infectious nodes. 

gsi{xs,xi,x R ) = 

^PkUk ^2 j x x s x : > I x R - l ~ : >Pr{d(S) = i,d{I) = j\ps,Pi}/ 

k 

J2P kUk J2 j xPr[d(S) = i,d{I)=j\p s ,pi] (14) 

k i,j\i+j<k 



d 

- — gs(xs,xi,x R ) 
axi 



-^—gs(xs,xi,x R ) 
dxi 



= g'(8(x s ps + xipi + x R (l - ps - Pi)))/g'{0) 



Although a rigorous proof for this is currently lacking, it is borne out by the 
success of this mathematical theory in predicting epidemic final size and dynamics 
(see sections |3j and 3.1 below). 
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Because arcs are distributed multinomially to nodes in sets S,I,1Z, we have 
gss(xs,xi,XR) = gsi(xs,xi,xn), which is easy to verify by repeating the cal- 
culation in equation |14| 

A useful property of PGF's is that the mean of the distribution they generate 
can be calculated by differentiating and evaluating with the dummy variables set 



to one [28]. Now using equations 13 and 14 we have the following results 

" d 



dx 



Ssi(I) = 
Ssi(S) = 



d 
dxi 

d 

dx s 



S31 = 
gsi(xs,xi,x R ) 

gsi(xs,xi,x R ) 



gsi(x,x,x) 



x s =x I =x R — l 



= eg"(9)/g'(9) (15) 

=i 

= Pi9g"(6)/g'{9) (16) 

= p s e g "(e)/g'(e) (17) 



J x s =x I —x R — l 

As a fraction — S nodes leave set S in time dt, the fraction of arcs between 
5 and I, Msi is reduced by the fraction of arcs from infectious nodes to the — S 
newly infectious nodes. Therefore M$i is reduced at rate — SSsi(I) /g' (1). Because 
Ssi(I) does not count the arc along which a node was infected, Msi is also reduced 
at a rate rM$i to account for all arcs which have an infectious ego which transmits 
to the susceptible alter. And in time dt, \il nodes become recovered. The average 
number of arcs in Ais per infectious node is proportional to Msi/I- Then Msi 
is reduced at a rate fiI{Msi / 1) = fiMsi- 

The quantity Msi is also increased, as new infected nodes have links to suscep- 
tible nodes. A newly infectious node will have on average Ssi{S) arcs to susceptible 
nodes, so Msi is increased at a rate — SSsi(S)/ g' (1). 

To summarize, Msi decreases at the sum of rates 

- -S6 S i(I)/g'(l) 

- rMsi 

- fJ-Msi 

And Msi increases at the sum of rates 



-SSsi(S)/g'(l) 
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Then applying equations |16| |17| and [TTJ we have 

Msi = {(-S)Ssi(S) - (~S)S S i(I))/g'(l) - (r + n)M S i 
= r Pl (p s - pi)6 2 g"(6)/g'(l) - (r + fj,)M S i 
Finally, it is necessary to determine the time derivative of Ms- 

Ms = j t e 9 '{e)/g'{l) = {Og' {6) + 66g"(6))/g'(l) 



(18) 



(19) 



= (-r Pl eg'(6)-r Pl e 2 g"(9))/g'(l) 

Now applying equations ^]9 18 and 19 to equation [8] we solve for pi in terms 
of the PGF and 6. 

pI = rpip s 0^Ml -pi(l-pi)r-pin (20) 

This equation makes use of the variable ps which changes in time. Deriving 
the dynamics of this variable will complete the model. This calculation is very 
similar to that for pj. 

. _ d M ss _ M ss M S M SS 

PS dt M s M s M 2 S [ 1 

The calculation for Mss is very similar to that for Msi- Newly infected nodes 
have on average Ssi{S) arcs to other susceptibles, so that 

M SS = -2x(-S)Ssi(S)/g'(l) 

(22) 

= -2rpip s 8 2 g"(e)/g'(l) 
where the factor of 2x accounts for two arcs per edge. 

Now applying equations [9] |19|and|22| to equation |21| we have 



Ps=rp lP s(l-e^l) (23) 



9'{0) 

The complete system of equations is summarized in table [3] 

3 The normalizing constant g'(l) cancels out and could have been left out these 
equations. 



SIR dyn. in rand. net. w. heter. conn. 



13 



Table 3. A summary of the nonlinear differential equations used to the describe 
the spread of a simple SIR type epidemic through a random network. The degree 
distribution of the network is generated by g(x). 



= -r Pl 6 

Pi = r PsPi°^^J ~ rPii 1 - Pi) - Pit* 

P* = r P3Pi j 1 - 9 ^) 

S = g(0) 

1 = r Pl 0g'(0)-»I 

The fraction of infectious nodes can be solved by introducing a fourth dynamic 
variable. The infectious class increases at a rate — S and decreases at a rate 
Therefore 

i = -r Pl 0g'(9)- f J,I (24) 

An advantage of dynamic modeling of epidemics in networks is that the time- 
evolution of variables besides incidence can be calculated. Above it was shown how 
to calculate the degree distribution among susceptible nodes (eqn. [T3|). Addition- 
ally, the degree distribution among nodes which are either infectious or recovered 
(set J~) can be calculated by taking the complement. 

gj(x) = (g(x)-g(ex))/(l-g(e)) (25) 
2.3. Initial Conditions 

If a small fraction e of the nodes in the network are selected uniformly at random 
and initially infected, we can anticipate the following initial conditions. The frac- 
tion of arcs with infectious ego will also be Mi = e, and since e is small, there is 
low chance of two initial infecteds being connected. Therefore Msi ~ Mi — e. 8, 
which can be interpreted as the fraction of degree one nodes remaining susceptible 
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will be 1 — e. And Ms = 1 — Msi — 1 — e because there are initially no recovered 
nodes. And Mss = Ms — Msi = 1 — 2e. To summarize, 

1. 9{t = 0) = 1 - e 

2. p/(t = 0) = Msi/Ms = e/(l - e) 

3. p s (t = 0) = Mss/Ms = (1 - 2e)/(l - e) 

2.^- Epidemic threshold 

Epidemic dynamics can fall into one of two qualitatively different regimes. Below 
a threshold in the ratio r/fj,, the final size (-Too) is necessarily proportional to the 
fraction of initial infectious nodes: loo oc e. But above this threshold, epidemics 
occur, and necessarily occupy a fraction of the population even as e — > 0. 

As per equation [1J the number of new infections in a small time interval is 
proportional to pi. This is in contrast to compartment models in which the number 
of new infections is proportional the current number of infectious. If pi(t — 0) < 0, 
an epidemic will necessarily die out without reaching a fraction of the population. 
The epidemic threshold occurs where 

pi(t = 0) = = TPsPidy^ ~ rp/(l - Pi) - PIU (26) 

Applying the initial conditions given in the last section and considering e < 1 
gives 

Pl(t = 0) = r \- 2 " 6 (1 - e)g"(9)/g'(6) - r^- 1 -^ - ji ' 



1-el-S < < l-el-e 

g"(0) \ 

y y ' r-fi ) = 



Rearranging yields the critical ratio r/fj, in terms of the PGF. 

w ">' - ^w^k (28) 

The epidemic threshold in equation [28] can also be put in terms of the the 
transmissibility, which is the probability that an infectious ego will transmit in- 
fection to a given alter. Integrating over an exponentially distributed duration of 
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infectiousness T, the transmissibility r is calculated to be 

f-OO 

t — Pr [transmit prior to T] x Pr [recover at T]dT 
Jt=o 

= f(l-e- tf )^ )dT=^- 
Jt=o r + ju 



(29) 



Then rearranging equation [28] yeilds the epidemic threshold in terms of r. 

r*= 5 '(l)/s"(l) (30) 
This is consistent with previous results based on bond-percolation theory |19] . 

3. Examples 

The model has been tested on several common degree distributions: 

k -z 

— Poisson: p k = ^-jn — . This is generated by 

g(x) = (31) 

— Power-law. For our experiments, we utilize power-laws with exponential cutoffs 
k: pt = , fc "! e _i/^, , k > 1 where Li n (x) is the nth polylogarithm of x. This is 
generated by 

g(x) = Li J (xe~ 1/K )/Li J (e~ 1/K ) (32) 

— Exponential: pk = (1 — e~ 1,/A )e~ Afc . This is generated by 

l-e" 1/A 

f W = 1 ZTfTX ( 33 ) 

1 — xe i 

Figure[l]shows the disease incidence for each of the degree distributions ( 31 1, ( 32 l 



and (|33|, with a force of infection r = .2 and recovery rate = .1. Initially 
e = 10 -4 nodes are infected. The parameters of the degree distributions were cho- 
sen so that each network has an identical average degree of 3. That is, the density 
of connections in each network is the same. Nevertheless, there is widely differ- 
ent epidemic behavior due to the different degree distributions. Consistent with 
previous research, the degree distribution has a great impact on the final size of 
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the epidemic |14II19] . More importantly, the three networks exhibit widely varying 
dynamical behavior. The power law network experiences epidemics which accel- 
erate very rapidly. Such epidemics enter the expansion phase (the time at which 
incidence increases at its maximum rate) virtually as soon as the first individual in 
the network is infected. Both the Poisson and exponential networks experience a 
lag before the expansion phase of the epidemic. These observations are consistent 
with the findings in [3] that the timescale of epidemics shortens with increasing 
contact heterogeneity. This has important implications for intervention strategies, 
as it is often the case that interventions are planned and implemented only after 
a pathogen has circulated in the population for some time. If an epidemic were to 
occur in the power law network, there would be little time to react before the the 
infection had reached a large proportion of the population. 

Recall from section [274| that below the epidemic threshold r* , only small, finite- 
sized outbreaks will occur. Figure [2] shows the qualitatively different dynamical 
behavior of outbreaks below the phase transition for networks with a Poisson 
distribution. Below the phase transition, the final size is always proportional to 
the fraction of initial infecteds e. 

Something offered by this model and not to the author's knowledge seen pre- 
viously, is an explicit calculation for how the degree distribution of susceptibles 
evolves over the course of the epidemic. We expect the degree distribution to 
become bottom-heavy, as high degree nodes are gradually weeded out of the pop- 
ulation of susceptibles. This is indeed observed in figure [3] for the Poisson trial 
described above. 



Recall that the degree distribution of susceptibles is generated by the multi- 
variate PGF (131. The explicit degree distribution can be retrieved from equa- 
tion [13] by differentiation. The following gives the probability that a susceptible 
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Fig. 1. The number of infecteds (including recovered) is shown versus time 
for an SIR model on three networks. Force of infection and mortality are con- 
stant: r = 0.2, fi — 0.1. The networks have Poisson (z — 3), power law 
(7 = 1.615, k = 20), and exponential (A = 3.475) degree distributions. Each 
of these degree distributions has an average degree of 3. 

node has m links at a time corresponding to 8. 




[~r-j;9s(x,x,x)] s =o/kl 



(34) 



For example, applying this to the Poisson PGF (equation (31 1 ) gives 



Pfe = 



(z6) k e 
fcT 



(35) 



which is simply the Poisson distribution with an adjusted parameter zxO. Another 
example is illustrated in figure [3] which shows the degree distribution among 
susceptibles for the power-law network considered above. 
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Fig. 2. The number of infecteds (including recovered) is shown versus time 
for an SIR model on a Poisson network {z = 3). Each of these trials are below 
the epidemic threshold required to sustain an epidemic. The outbreak size is re- 
ported as a multiple of the fraction of initial infecteds in the network. Mortality 
is constant, /i = 0.4, while three different levels of the force of infection are tried, 
r = 0.15,0.17,0.18. 

3.1. Stochastic Simulations 

Simulation of SIR on networks presents two challenges: A random network must 
be generated with the desired degree distribution. Secondly, the stochastic rules 
that govern the transmission of disease at the microscopic scale must be well- 
defined, and an algorithm must be developed to aggregate this behavior into a 
large-scale simulation. 

The random generation of networks with a given degree distribution is a well- 
explored problem. The first algorithm was proposed by Molloy and Reed [17] 
which I have used for these experiments. Subsequent research has shown that 
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Degiee 



Fig. 3. The degree distribution for susceptible nodes where the epidemic size 
is 50%, 75%, and 100% of the final size, as well as degree distribution at the 
beginning of the epidemic. The degree distribution for the network as a whole is 



a power law with exponential cutoff (equation 32 1 



imperfections can arise in the networks generated by this algorithm, but such 
biases should be tolerably small for these purposes [16] . 
The simulation dynamics are as follows: 

— A node is chosen uniformly at random from the network as an initial infected. 

— An infected node v will recover after an exponentially distributed random time 
interval At^ ~ Exp(ji). 

— When a node v is infected, each arc (v, x) has a time of infection At x drawn 
from an exponential distribution Exp(r). If At x < At^, node x is infected 
after time At x . Otherwise x is not infected by v. 

This process continues until there are no more infectious nodes. 

Figure [4] shows the results of 450 simulations for the Poisson random network 
considered in the last section (z = 3) with 10 4 nodes. The black dotted line 
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9000 




10 20 30 40 50 60 70 80 90 100 

Time (t) 

Fig. 4. 450 simulation trajectories of the cumulative epidemic incidence J (dot- 
ted lines) for a Poisson (z = 3) random network. The solid blue line shows the 
analytical solution. 

represents an independent simulation trajectory. The thick, blue line that cuts 
through the dense mass of simulation trajectories is the analytical trajectory based 
on the equations in table [3] The initial conditions were chosen as in the previous 
section using e = 10~ 4 . 

Figure [5] shows a similar series of simulations for the power law degree dis- 
tribution considered in the last section. In both cases, the analytical trajectory 
traverses the region with the highest density of simulation trajectories. The simu- 
lation trajectories also exhibit significant variability in the time required to reach 
the expansion phase and final size. This is largely due to the significant impact 
of random events early on in the epidemic. For example, an initial infected with 
a low average degree, or one which takes an inordinate amount of time to infect 
the next infected can markedly delay the onset of the expansion phase. 

Figure [6] shows the median-time incidence for the exponential and Poisson 
networks discussed in the last section. The data points show the median time 
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Fig. 5. 450 simulation trajectories of the cumulative epidemic incidence J for a 
power law (7 = 1.615, k = 20) random network. The solid line shows the analytical 
solution based on the system of equations in table [3] 

required to reach a given incidence among 450 simulation trajectories. The solid 
line shows the analytical trajectory based on the system of equations given in 
table [3] Intuitively, the data points are showing the path of the most central 
trajectory from the swarm of simulation trajectories such as in figure [4] 



4. Discussion 



The statistical properties of SIR epidemics in random networks have been un- 
derstood for some time, but the explicit dynamics have been understood mainly 
through simulation. This paper has addressed this shortcoming by proposing a 
system of nonlinear ordinary differential equations to model SIR dynamics in ran- 
dom networks. 
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A Exponential Median Time Incidence 
Pojsson Median Time Incidence 



1,0 



Time (t) 



Fig. 6. The median time required to reach a given incidence J is shown for a 
Poisson network (z = 3, circles) and an exponential network (A = 3.475, triangles). 
The solid line shows the analytical solution based on the system of equations in 
table [U 

It should be noted that the SI dynamics are a special case of this model 
(fi = 0), in which case the ultimate extent of the epidemic is simply the giant 
component of the network. 

The distribution of contacts, even holding the density of contacts constant, 
has enormous impact on epidemic behavior. This goes beyond merely the extent 
of the epidemic, but as shown here, the dynamical behavior of the epidemic. In 
particular, the distribution of contacts plays a key role in determining the onset 
of the expansion phase. 

The distribution dynamics from equation [13] and shown in figure [3] have im- 
portant implications for vaccination strategies. Previous work 12, 10] has focused 

4 The giant component of a network, if it exists, is a set of nodes such there 
exists a path between any two of nodes, and furthermore occupies a non-zero 
fraction of the network in the limit as network size goes to infinity. 
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on determining the critical levels of vaccination required to halt or prevent an 
epidemic. It is usually taken for granted that contact patterns among suscepti- 
bles are constant. Furthermore, most widespread vaccinations occur only once an 
epidemic is underway. Future research could be enhanced by considering optimal 
vaccination levels when the epidemic proceeds unhindered for variable amounts 
of time. 

It is hoped that the distribution dynamics described in this paper will find 
applications beyond modeling heterogeneous connectivity. The dynamic PGF ap- 
proach might be used to capture other forms of heterogeneity, such as of suscep- 
tibility, mortality, and infectiousness. 
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